#include "AnaDef.h"
#include "style.h"

const Double_t EPSILON = 1e-10;

void stylec0()
{
  style::fgkTitleSize = 0.07;
  style::fgkTextSize = 0.07;//for legend
  style::fgkYTitleOffset = 1;
}

void stylec1()
{
  style::fgkTitleSize = 0.1;//title
  style::fgkTextSize = 0.1;//label, legend (if any)
  style::fgkYTitleOffset = 0.7;
  style::fgkXTitleOffset = 1.35;
}

//3000
Double_t IDdEdxM3Q0(const Double_t *xx, const Double_t *par)
{
  const Double_t mean=xx[0];
  return                    TMath::Abs(par[0])/TMath::Sqrt(mean) + par[1] + par[2]/mean;
}

//2000
Double_t IDdEdxM2Q0(const Double_t *xx, const Double_t *par)
{
 const Double_t mean=xx[0];
  return                  ( TMath::Abs(par[0]) * TMath::Power(mean,par[1]) + par[2] )/mean;
}

Double_t RawAleph(const Double_t *xx, const Double_t *par)
{
  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = TMath::Abs(par[0]);
  const Double_t p1 = TMath::Abs(par[1]);
  const Double_t p2 = par[2];//TMath::Abs(par[2]);//still first have to be positive
  const Double_t p3 = TMath::Abs(par[3]);
  const Double_t p4 = TMath::Abs(par[4]);//has to be positive so that in bg-> inf there is no bg dependence anymore

  //return p0*((p1-TMath::Log(p2+1/TMath::Power(bg,p4)))/TMath::Power(beta,p3) - 1 );
  //return p0*( (p1- TMath::Log(p2+EPSILON)- TMath::Log(1+1/(TMath::Power(bg,p4)*p2+EPSILON))  )/( TMath::Power(beta,p3)+EPSILON ) - 1 );
  //log(p2) -> p2 because p2 sometimes ~EPSILON
  return p0*( (p1- p2 - TMath::Log(1+1/(TMath::Power(bg,p4)*TMath::Exp(p2)+EPSILON))  )/( TMath::Power(beta,p3)+EPSILON ) - 1 );
  //p2 can be very small ~ 1e-9, do not separate p2 out! -> absorb p2 into p1, remember to change initial value of p1
  //return p0*( (p1- TMath::Log(1+1/(TMath::Power(bg,p4)*p2+EPSILON))  )/( TMath::Power(beta,p3)+EPSILON ) - 1 );
  //change p0 as plateau
  //return p0/(p1-1)*( (p1- TMath::Log(1+1/(TMath::Power(bg,p4)*p2+EPSILON))  )/( TMath::Power(beta,p3)+EPSILON ) - 1 );
  //redefine p1 so that p1>1 --> p1>0
  //return p0/p1*( (p1 + 1 - TMath::Log(1+1/(TMath::Power(bg,p4)*p2+EPSILON))  )/( TMath::Power(beta,p3)+EPSILON ) - 1 );
}

Double_t myBetheBloch(const Double_t *xx, const Double_t *par)
{
  //RawALEPH has numerical difficulty in ln(p2+1/pow(bg, p4))
  //use flat dedx for electron. deviation from this at high bg is not inportant because anyway there the statistics is very limited
  //npar = 5;

  const Double_t bg = xx[0];
  const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  const Double_t p0 = par[0];//overall prefactor
  const Double_t p1 = par[1];//beta exponent ~ 2
  const Double_t p2 = par[2];
  const Double_t p3 = par[3];
  const Double_t p4 = par[4];//plateau ~ 80

  //EPSILON added to protect 0 denominator!
  const Double_t bethebloch = p0 *( 1 + p2*TMath::Log(bg) + p3*TMath::Log(beta) )/( TMath::Power(beta, p1)+EPSILON) ;

  //tp = 100, electron momentum = 51 MeV, pion momentum = 14 GeV
  const Double_t tp = 10;

  Double_t qq = 0;
  if(bg<tp){
    qq = bethebloch;
  }
  else{
    qq = TMath::Min(p4, bethebloch);
  }

  return qq;
}

TF1 * GetmyBetheBloch(const TString tag, const Double_t pars[], const Int_t col, const Int_t sty, const Double_t xmin, const Double_t xmax)
{
  TF1 *f1=new TF1(Form("myBetheBloch%s", tag.Data()), myBetheBloch, xmin, xmax, 5);
  f1->SetLineColor(col);
  f1->SetLineStyle(sty);
  f1->SetParameters(pars);
  return f1;
}

TF1 * GetRawAleph(const TString tag, const Double_t pars[], const Int_t col, const Int_t sty, const Double_t xmin, const Double_t xmax)
{
  TF1 *f1=new TF1(Form("RawAleph%s", tag.Data()), RawAleph, xmin, xmax, 5);
  f1->SetLineColor(col);
  f1->SetLineStyle(sty);
  f1->SetParameters(pars);
  return f1;
}

TF1 * GetRes3000(const TString tag, const Double_t pars[], const Int_t col, const Int_t sty, const Double_t xmin, const Double_t xmax)
{
  TF1 *f1=new TF1(Form("Res3000%s", tag.Data()), IDdEdxM3Q0, xmin, xmax, 3);
  f1->SetLineColor(col);
  f1->SetLineStyle(sty);
  f1->SetParameters(pars);
  return f1;
}

TF1 * GetRes2000(const TString tag, const Double_t pars[], const Int_t col, const Int_t sty, const Double_t xmin, const Double_t xmax)
{
  TF1 *f1=new TF1(Form("Res2000%s", tag.Data()), IDdEdxM2Q0, xmin, xmax, 3);
  f1->SetLineColor(col);
  f1->SetLineStyle(sty);
  f1->SetParameters(pars);
  return f1;
}

void drawFunc(TPad * pad0, TPad * pad1, TF1 *ff, TF1 *fref, const TString dopt, const TString ytit, const Double_t ymin, const Double_t ymax, const Double_t rmin, const Double_t rmax)
{
  stylec0();
  style::PadSetup(pad0);
  pad0->SetLeftMargin(0.14);
  pad0->SetRightMargin(0.03);
  pad0->SetTopMargin(0.1);
  pad0->SetBottomMargin(0.01);
  pad0->cd();

  //histogram style works normal
  TH1D * h0 = style::Func2Hist(ff, kTRUE);
  style::ResetStyle(h0);
  h0->SetMinimum(ymin);
  h0->SetMaximum(ymax);
  h0->SetLineWidth(2);
  h0->Draw(dopt);

 //----------------------
  stylec1();
  style::PadSetup(pad1);
  pad1->SetTopMargin(0.01);
  pad1->SetLeftMargin(0.14);
  pad1->SetRightMargin(0.03);
  pad1->SetBottomMargin(0.3);
  pad1->cd();

  //GetHistogram very very bad
  TH1D * hf = style::Func2Hist(ff, kTRUE); //(TH1D*)ff->GetHistogram()->Clone(Form("%sff",ff->GetName()));
  TH1D * href = style::Func2Hist(fref, kTRUE);//fref->GetHistogram();

  hf->Divide(href);
  style::ResetStyle(hf);

  hf->SetMinimum(rmin);
  hf->SetMaximum(rmax);
  hf->GetXaxis()->SetTitle(ff->GetXaxis()->GetTitle());
  hf->GetYaxis()->SetTitle(ytit);
  hf->GetYaxis()->SetNdivisions(505);
  hf->SetLineWidth(2);
  hf->Draw(dopt);

}

void drawdEdxCompJetInc(const Int_t kMC=0)
{
  style::SetGlobalStyle(0);

  TCanvas *cc=new TCanvas("cc","",600,500);

  TPad *c0=new TPad("c0","",0,0.4,1,1);
  c0->SetLogx(1);
  c0->Draw();

  TPad *c1= new TPad("c1","", 0,0,1,0.4);
  c1->SetLogx(1);  
  c1->Draw();

  //-----------------------------------
  Double_t fPars[10];

  //DATA
  //xlu@lxsub01:/hera/alice/xlu/workspace/hera2dbetabloch/FracFunc_10-d-e/300by2400_3tpc_jet_12122013/finish$ for ii in $(echo TPC1TOF0/*ETA1*_{0-100,15-20}*/XB5YB6*aleph*anchor4*/ID3*/logmn2d.txt ); do echo $ii ; cat $ii | grep -vi bin | tail -n 20 | head -n 11;done

  //MC
  //xlu@lxsub01:/hera/alice/xlu/workspace/hera2dbetabloch/FracFunc_10-d-e/MC_300by2400_20131208/finish$ for ii in $(echo TPC1TOF0/*ETA1*_{0-100,15-20}*/XB5YB6*aleph*anchor4*/ID3*/logmn2d.txt ); do echo $ii ; cat $ii | grep -vi bin | tail -n 20 | head -n 11;done

  if(kMC==0){
    /*
TPC1TOF0/DATA.histRecTHnIncl_LHC10de_TPC1TOF0_kMC0_ETA1_JPTall_0-100.coherentr04/XB5YB6_aleph_anchor4_RegOpt0/ID3000_2_NTYPE4FixHighPElectron/logmn2d.txt
fNdata = 14400
  Status: fkErr = 0, fStatChi2 = 1.058365e+04, NDF() = 14254, ReducedChi2() = 7.425037e-01, fBeta =  1.000000e+00, fPenalty = 1.531539e+01, b*p = 1.531539e+01 fcn = 1.059896e+04
    */
    fPars[0] = 1.329709e+00;
    fPars[1] = 2.943341e+01;
    fPars[2] = -2.818210e+01;
    fPars[3] = 2.066739e+00;
    fPars[4] = 4.859660e+00;
    fPars[5] = 3.814119e-07;
    fPars[6] = 5.257955e-02;
    fPars[7] = 1.137848e-01;
  }
  else{
    /*
TPC1TOF0/DATA.histRecTHnIncl_LHCMC10de_TPC1TOF0_kMC1_ETA1_JPTall_0-100.coherentr04/XB5YB6_aleph_anchor4_RegOpt0/ID3000_2_NTYPE4FixHighPElectron/logmn2d.txt
fNdata = 14000
  Status: fkErr = 0, fStatChi2 = 9.505574e+03, NDF() = 13857, ReducedChi2() = 6.859763e-01, fBeta =  1.000000e+00, fPenalty = 1.510291e+01, b*p = 1.510291e+01 fcn = 9.520677e+03
    */
    fPars[0] = 1.469075e+00;
    fPars[1] = 2.412253e+01;
    fPars[2] = -2.732735e+01;
    fPars[3] = 1.953326e+00;
    fPars[4] = 4.606017e+00;
    fPars[5] = -8.931372e-07;
    fPars[6] = 4.659478e-02;
    fPars[7] = 3.239762e-01;
  }

  TF1 * fmeanInc=GetRawAleph("meanInc",fPars, kBlack, kDashed, 0.4, 1e4);
  fmeanInc->GetXaxis()->SetTitle("#beta#gamma");
  fmeanInc->GetYaxis()->SetTitle("#mu (arb. units)");

  TF1 * fsigmaInc=GetRes3000("sigmaInc", &(fPars[5]), kBlack, kDashed, 40, 90);
  fsigmaInc->GetXaxis()->SetTitle("#mu (arb. units)");
  fsigmaInc->GetYaxis()->SetTitle("#sigma'");


  //-----------------

  if(kMC==0){
    /*
TPC1TOF0/DATA.histRecTHnIDFF_LHC10de_TPC1TOF0_kMC0_ETA1_JPT4_15-20.coherentr04/XB5YB6_aleph_anchor4_RegOpt0/ID3000_2_NTYPE4FixHighPElectron/logmn2d.txt
fNdata = 12000
  Status: fkErr = 0, fStatChi2 = 6.295890e+03, NDF() = 11872, ReducedChi2() = 5.303142e-01, fBeta =  1.000000e+00, fPenalty = 3.276327e+00, b*p = 3.276327e+00 fcn = 6.299166e+03
    */
    fPars[0] = 1.359865e+00;
    fPars[1] = 2.886111e+01;
    fPars[2] = -2.766997e+01;
    fPars[3] = 2.093354e+00;
    fPars[4] = 4.751080e+00;
    fPars[5] = -4.964920e-06;
    fPars[6] = 5.979800e-02;
    fPars[7] = -3.438665e-01;
  }
  else{
    /*
TPC1TOF0/DATA.histRecTHnIDFF_LHCMC10de_TPC1TOF0_kMC1_ETA1_JPT2_15-20.coherentr04/XB5YB6_aleph_anchor4_RegOpt0/ID3000_2_NTYPE4FixHighPElectron/logmn2d.txt
fNdata = 12000
  Status: fkErr = 0, fStatChi2 = 5.404027e+03, NDF() = 11872, ReducedChi2() = 4.551910e-01, fBeta =  1.000000e+00, fPenalty = 5.105324e+00, b*p = 5.105324e+00 fcn = 5.409132e+03
    */
    fPars[0] = 1.475429e+00;
    fPars[1] = 2.420627e+01;
    fPars[2] = -2.704660e+01;
    fPars[3] = -1.955279e+00;
    fPars[4] = 4.541755e+00;
    fPars[5] = 2.828524e-07;
    fPars[6] = 5.552650e-02;
    fPars[7] = -1.636980e-01;
  }

  TF1 * fmeanJet=GetRawAleph("meanJet", fPars, kRed, kSolid, 0.4, 1e4);
  fmeanJet->GetXaxis()->SetTitle("#beta#gamma");
  fmeanJet->GetYaxis()->SetTitle("#mu (arb. units)");

  TF1 * fsigmaJet=GetRes3000("sigmaJet", &(fPars[5]), kRed, kSolid, 40, 90);
  fsigmaJet->GetXaxis()->SetTitle("#mu (arb. units)");
  fsigmaJet->GetYaxis()->SetTitle("#sigma'");
  
  //========================================================

  TLegend *lg=new TLegend(0.55, 0.65, 0.98, 0.88);
  stylec0();
  style::ResetStyle(lg);

  drawFunc(c0,c1, fmeanInc, fmeanInc, "l", "ratio to ref.", 45, 135, 0.995, 1.005);
  drawFunc(c0,c1, fmeanJet, fmeanInc, "l same", "",         45, 135, 0.995, 1.005);

  c0->cd();

  lg->SetHeader(Form("%s10d+e 0.2<|#eta|<0.4", kMC?" MC":""));
  lg->AddEntry(fmeanInc,"inclusive (ref.)","l");
  lg->AddEntry(fmeanJet,"jets pt 15-20 GeV/c", "l");
  lg->Draw();

  cc->Print(Form("meanIncJetkMC%d.eps", kMC));
  //========================================================

  lg=new TLegend(0.55, 0.65, 0.98, 0.88);
  stylec0();
  style::ResetStyle(lg);

  drawFunc(c0,c1, fsigmaInc, fsigmaInc, "l", "ratio to ref.",        0.039, 0.081, 0.88, 1.12);
  drawFunc(c0,c1, fsigmaJet, fsigmaInc, "l same", "",     0.039, 0.081, 0.88, 1.12);

 
  //c0->SetTopMargin(0.08);
  c0->SetLogx(0);
  c1->SetLogx(0);
  c0->cd();

  lg->SetHeader(Form("%s10d+e 0.2<|#eta|<0.4", kMC?" MC":""));
  lg->AddEntry(fsigmaInc,"inclusive (ref.)","l");
  lg->AddEntry(fsigmaJet,"jets pt 15-20 GeV/c", "l");
  lg->Draw();

  cc->Print(Form("sigmaIncJetkMC%d.eps", kMC));
}

int main()
{
  drawdEdxCompJetInc();
  return 0;
}
